# Timing Matters: Climate Policies and Adaptive Resilience
# Visualizations
# D'Amico & Maboudi (2025)

# ============================================================================
# SETUP
# ============================================================================

rm(list = ls())

library(changepoint)
library(dplyr)
library(ggplot2)
library(patchwork)
library(readr)
library(rnaturalearth)
library(sf)
library(showtext)

# Load cleaned data
pdata <- read_csv("pdata_clean.csv")

# Load RF results
rf_results <- readRDS("rf_variable_importance.rds")

# Set up Times New Roman font (adjust path as needed)
font_add(family = "Times New Roman", 
         regular = "/Library/Fonts/Times New Roman.ttf")
showtext_auto()

# ============================================================================
# FIGURE A.1: CHANGE POINT ANALYSIS (APPENDIX A)
# ============================================================================

# Aggregate climate law features by year
climate_features <- pdata %>%
  filter(!is.na(year)) %>%
  group_by(year) %>%
  summarise(
    across(c(adaptation_executive_count, adaptation_legislative_count,
             adaptation_plans_count, cl_domestic_count), 
           ~ sum(.x, na.rm = TRUE)),
    .groups = 'drop'
  ) %>%
  mutate(total_features = adaptation_executive_count + adaptation_legislative_count + 
           adaptation_plans_count + cl_domestic_count)

# Calculate aggregated percent change
climate_features <- climate_features %>%
  arrange(year) %>%
  mutate(agg_pct_change = (total_features - lag(total_features)) / 
           (lag(total_features) + 1) * 100)

# Change point detection
cpt_data <- na.omit(climate_features$agg_pct_change)
cpt_model <- cpt.mean(cpt_data, method = "PELT")

# Plot
fig_a1 <- ggplot(climate_features, aes(x = year, y = agg_pct_change)) +
  geom_line(size = 1) +
  geom_vline(xintercept = c(1998, 2008, 2016), 
             linetype = "dashed", color = "red") +
  labs(x = "Year", 
       y = "Aggregated % Change in Climate Law Features") +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"))

ggsave("Figure_A1_changepoint.png", fig_a1, width = 8, height = 5, dpi = 300)

# ============================================================================
# FIGURE A.2: RANDOM FOREST VARIABLE IMPORTANCE (APPENDIX E)
# ============================================================================

# Create variable labels
variable_labels <- c(
  "adaptation_law_stock" = "Adaptation Law Stock",
  "imports_pGDP" = "Imports (% of GDP)",
  "max_populism" = "Maximum Populism",
  "ROL_perc" = "Rule of Law (%)",
  "percent_gdp_services" = "Services (% of GDP)",
  "auton" = "Federalism Index",
  "disaster_count" = "Disaster Count",
  "ln_gdp" = "Log GDP",
  "gdpgrowth_per" = "GDP Growth Rate (%)",
  "adaptationaid_recipient" = "Adaptation Aid Recipient",
  "envexp_pgdp" = "Environmental Expenditure (% GDP)",
  "global_north" = "Global North",
  "financial_recession" = "Financial Recession",
  "adaptation_plans_stock" = "Adaptation Plans (Stock)",
  "ln_gdp_square" = "Log GDP Squared"
)

# Function to create RF plot
create_rf_plot <- function(rf_result) {
  imp_data <- rf_result$importance_data %>%
    slice_head(n = 15) %>%
    mutate(Variable_Label = ifelse(Variable %in% names(variable_labels),
                                    variable_labels[Variable],
                                    Variable))
  
  ggplot(imp_data, aes(x = reorder(Variable_Label, `%IncMSE`), y = `%IncMSE`)) +
    geom_col(fill = "#CBC3E3") +
    coord_flip() +
    labs(subtitle = paste("Temporal Span:", rf_result$segment),
         x = NULL,
         y = "% Increase in MSE") +
    theme_minimal() +
    theme(
      plot.subtitle = element_text(hjust = 0.5, size = 10, family = "Times New Roman"),
      axis.text = element_text(size = 8, family = "Times New Roman"),
      axis.title.x = element_text(size = 10, family = "Times New Roman"),
      text = element_text(family = "Times New Roman")
    )
}

# Create plots for each temporal segment
rf_plots <- lapply(rf_results, create_rf_plot)

# Combine into single figure
fig_a2 <- wrap_plots(rf_plots, ncol = 1)

ggsave("Figure_A2_rf_importance.png", fig_a2, width = 8, height = 12, dpi = 300)

# ============================================================================
# FIGURE A.3: TRENDS IN LAW STOCK AND CAPACITY (APPENDIX I)
# ============================================================================

# Convert year to numeric
pdata$year <- as.numeric(as.character(pdata$year))

# Filter from 1995 onwards
pdata_filtered <- pdata %>% filter(year >= 1995)

# Aggregate by year
law_agg <- pdata_filtered %>%
  group_by(year) %>%
  summarise(
    mean_law = mean(adaptation_law_stock, na.rm = TRUE),
    se_law = sd(adaptation_law_stock, na.rm = TRUE) / sqrt(sum(!is.na(adaptation_law_stock))),
    .groups = 'drop'
  )

capacity_agg <- pdata_filtered %>%
  group_by(year) %>%
  summarise(
    mean_capacity = mean(capacity, na.rm = TRUE),
    se_capacity = sd(capacity, na.rm = TRUE) / sqrt(sum(!is.na(capacity))),
    .groups = 'drop'
  )

# Plot 1: Law stock
plot_law <- ggplot(law_agg, aes(x = year, y = mean_law)) +
  geom_line(color = "purple", size = 1) +
  geom_errorbar(aes(ymin = mean_law - se_law, ymax = mean_law + se_law),
                color = "purple", width = 0.3, size = 0.5) +
  labs(x = "Year", y = "a. Adaptation Law Stock") +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"))

# Plot 2: Capacity
plot_capacity <- ggplot(capacity_agg, aes(x = year, y = mean_capacity)) +
  geom_line(color = "red", size = 1) +
  geom_errorbar(aes(ymin = mean_capacity - se_capacity, 
                    ymax = mean_capacity + se_capacity),
                color = "red", width = 0.3, size = 0.5) +
  labs(x = "Year", y = "b. Average Capacity") +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"))

# Combine
fig_a3 <- plot_law + plot_capacity + plot_layout(ncol = 2)

ggsave("Figure_A3_trends.png", fig_a3, width = 10, height = 4, dpi = 300)

# ============================================================================
# FIGURE A.4: POPULISM TRENDS (APPENDIX J)
# ============================================================================

populism_agg <- pdata %>%
  filter(year >= 1960, year <= 2023) %>%
  group_by(year) %>%
  summarise(avg_populism = mean(max_populism, na.rm = TRUE), .groups = 'drop')

fig_a4 <- ggplot(populism_agg, aes(x = year, y = avg_populism)) +
  geom_line(color = "purple", size = 1) +
  labs(x = "Year", 
       y = "Average Populism",
       title = "Average Levels of Populism Over Time (1960-2023)") +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"))

ggsave("Figure_A4_populism.png", fig_a4, width = 8, height = 5, dpi = 300)

# ============================================================================
# FIGURE A.5: CLIMATE EVENTS OVER TIME (APPENDIX K)
# ============================================================================

climate_events <- pdata %>%
  filter(year >= 1992, year <= 2021) %>%
  group_by(year) %>%
  summarise(total_events = sum(exp(disaster_count) - 1e-6, na.rm = TRUE), 
            .groups = 'drop')

fig_a5 <- ggplot(climate_events, aes(x = year, y = total_events)) +
  geom_line(color = "purple", size = 1) +
  labs(x = "Year",
       y = "Total Climate Events",
       title = "Total Climate Events Over Time (1992-2021)") +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"))

ggsave("Figure_A5_climate_events.png", fig_a5, width = 8, height = 5, dpi = 300)

# ============================================================================
# MAPS: LAW STOCK AND CAPACITY (2023)
# ============================================================================

# Load world map
world <- ne_countries(scale = "medium", returnclass = "sf")

# Filter for 2023
pdata_2023 <- pdata %>% filter(year == 2023)

# Merge with map data
merged_data <- world %>%
  left_join(pdata_2023, by = c("iso_a3" = "iso3c"))

# Map 1: Adaptation law stock
map_laws <- ggplot() +
  geom_sf(data = merged_data, aes(fill = adaptation_law_stock)) +
  scale_fill_gradient2(name = "a. Adaptation Law Stock", 
                       low = "blue", mid = "white", high = "red", 
                       midpoint = median(merged_data$adaptation_law_stock, na.rm = TRUE)) +
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(family = "Times New Roman"))

# Map 2: Capacity
map_capacity <- ggplot() +
  geom_sf(data = merged_data, aes(fill = capacity)) +
  scale_fill_gradient2(name = "b. Capacity Score", 
                       low = "blue", mid = "white", high = "red", 
                       midpoint = median(merged_data$capacity, na.rm = TRUE)) +
  theme_void() +
  theme(legend.position = "bottom",
        text = element_text(family = "Times New Roman"))

# Combine maps
maps_combined <- map_laws + map_capacity + plot_layout(nrow = 2)

ggsave("Figure_maps.png", maps_combined, width = 10, height = 8, dpi = 300)

# ============================================================================
# SUMMARY
# ============================================================================

cat("\n=== VISUALIZATIONS COMPLETE ===\n")
cat("Figures saved:\n")
cat("  - Figure_A1_changepoint.png\n")
cat("  - Figure_A2_rf_importance.png\n")
cat("  - Figure_A3_trends.png\n")
cat("  - Figure_A4_populism.png\n")
cat("  - Figure_A5_climate_events.png\n")
cat("  - Figure_maps.png\n")
